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ABSTRACT 

This  research  contribution  describes  a  computer  program  (CNA  Number 
76-67)  which  determines  the  Discrete  Fourier  Transform  of  a  set  of  data,  using 
a  recently  developed  technique  known  as  the  Fast  Fourier  Transform.  The  re¬ 
lation  between  Discrete  Fourier  Transforms  and  Fourier  Series  when  the  data 
is  periodic  is  also  shown. 
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in  the  derivations),  which  will  usually  entail  some  extra  input  and  output  operations 
to  keep  the  format  in  line  with  the  format  in  the  derivations. 

Appendix  A  contains  the  definition  of  the  Discrete  Fourier  Transform  and 
its  relation  to  the  Fourier  Series.  Appendix  B  is  a  description  of  the  computer 
program,  which  uses  the  FFT  technique,  along  with  its  limitations  and  some 
possible  uses.  Also  included  is  a  listing  of  the  program  itself,  which  is  written 
in  FORTRAN  II  for  the  CDC  3400  computer.  Appendix  CJ  gives  an  example  of  how 
the  program  can  be  used  and  some  numerical  results. 
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APPENDIX  A 


APPENDIX  A 


DEFINITIONS 


i liven  a  time  series  x(ka  ).  k=H . 

v\\\  rsc  arc  tUfiiu  d  as  follr.w 

N-l  ,  the  Discrete  Fourier  Transform  (I\ 

’  F)  and  its 

N-l 

51  x(k)exp  (  -7r‘fkn) 
k-0 

n=0,  1 . N-l 

(  I  l) 

N-l 

and  x(k)=  5Z  H(n)exp  ). 

n-o 

k=(),  1 . N-l 

(lb) 

where  i 

=  Vm 

Bin)  will  in  general  bi  complex  when  x  is  real.  The  similarity  lie  tween  the  DFT  and  the  Fourier 
Series  is  evident  when  we  consider  the  exponential  form  of  the  Fourier  Series: 

C<n)=i  J  x(t)cxp  (  )  dt, 

1  0 

all  integer  n. 

(2a) 

n=* 

and  x(t>  51  C(n)exp  ( -l~~  j‘nt )  . 

n=-a. 

where  x(t)  is  piecewise  continuous  and 
periodic  with  period  T. 

(2b) 

Now  if  the  above  integral  is  approximated  by  its  Riemann  sum,  C(n)  becomes  approximately 

S’- 1 

C(n)  ^  5-  x(kAt>cxp  (  where  T=Nat.  In  other  words  we  use  N 

k=0  ‘  samples  in  the  semi-open  interval  [0.  T). 

Hut  th  sum  on  the  right  is  just  R(n).  so  we  liavc 

C'(n)=B(n>  ( 3) 

I  Ills  then  is  tlie  relation  between  Fourier  Scries  and  Discrete  Fourier  Transforms  when  x(t)  is 
periodic.  Tlie  (  (n)  tell  us  how  much  of  x(t)  can  be  attributed  to  sinusoids  of  "frequency" 

,=T=K  (if)' 

and  can  be  approximated  by  the  DFT  coefficients  B(n).  This  becomes  clearer  when  we  realize 
that  in  fact  (reference  (hi) 

<r 

li(n)=52  l  <n+|.\)  .  n=0,  1 . N-l.  (4) 

J”-® 

I  hit  is.  IK  it)  is  the  sum  of  overlapped  segments  of  C(n>.  Figure  A-l  shows  this  relationship 
N  tween  Bin)  and  (’till.  In  order  to  make  il(n)  a  better  approximation  to  C(n)  we  must  increase 
ilu-  number  of  samples  in  the  pi  riod  T. 

Now  if  x(k)  is  real,  C(n)  will  be  an  even  function  of  n;  that  is,  C(n)=T(-n).  F.quation  (4)  thi  n 
giws  ns  some  further  information  about  B(n) 
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FIG.  A— 1:  SHOWING  HOW  THE  OFT  COEFFICIENTS  ARE  RELATED  TO  THE 
FOURIER  SERIES  COEFFICIENTS  OF  A  TIME  SERIES 
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co  x 

B(N-n)=  £  C(  -n+jN+N)  =  22  C(-n-Kj+l)N) 
j=-»  j=-oc 

00 

-  £  C(-n-t-mN)  with  m=j+l 
m=-» 

00 

=  2-  C(n-mN)  since  C  is  even 

m=-a> 

-00 

=  22  C(n+kN)  with  k=-m 

k=oo 

to 

=  22  <■  (n+kN)  since  the  order  of  summation  is  immaterial 

k=-oo 

B(N-n)=B(n)  n=0,  1 . N-l 

Thus  B(n)  is  symmetrical  about  n=N/2  (see  figure  A-l).  Care  must  be  taken,  therefore, 
not  to  use  B(n)  as  an  approximation  to  C(n)  when  hS-N/2. 


(5) 
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APPENDIX  B 


DESCRIPTION  OF  TMF  PROCR  AM 

Tin-  con  of  iik  program  to  copipur.  the  DFT  of  a  series  is  really  quite  compact,  and  ..  j, 
tit  express, si  a  a  ’o’  i  nesteii  do  loops: 

L  )  ]  mi  j  ] ,  s 
DO  100  k"  2:  ~J-1 
DO  100  n=0,  2J_I-1 
IXilOOnnO  i 

lot*  BlnHc- 2^ h;t  2^  *)-H(n4k-  2-*)-H-l)mB(n-He  2^+2^  *)exp 

where  N=2S  data  jxiints  are  used.  * 

An  advantage  of  the  program  not  mentioned  in  the  body  of  the  text  is  that  the  c<  efficients  can  use  the 
same  storage  as  the  original  data,  since  the  program  exchanges  pairs  of  values  (m=0,  1  in  the  above 
formulation)  after  appropriately  weighting  them.  An  important  programming  consideration,  however 
is  that  all  arrays  must  he  in  complex  notations  to  effect  this  space-saving,  so  real  time  series  must 
he  converted  to  a  complex  format  before  being  used. 

The  program  allows  up  to  1024=2^  data  points.  To  store  up  to  2*^  points,  the  user  need 
only  redimension  the  arrays  to  set  aside  that  much  storage.  If  still  more  data  is  used,  then 
i  (implicated  changes  must  he  made  to  the  subroutine  written  in  COMPASS. 

The  user  may  also  call  for  the  inverse  transform  to  get  back  a  time  series  from  a  set  of  co- 
,  llicients  by  calling  the  subroutine  INVERSE.  At  any  time  after  ceiling  the  subroutines  for  tin 
tr  insform  or  inverse  transform,  the  results  arc  stored  (in  complex  form)  in  the  original  data 
locations,  available  to  the  user  for  printout  or  manipulation. 

As  a  final  feature,  the  user  may  call  subroutines  10  smooth  the  coefficients.  Thu  reason  for 
wanting  to  smooth  the  DFT  coefficients  is  that  our  data  extends  only  over  a  finite  time  interval; 
this  however,  usually  only  represents  the  portion  of  the  process  we  have  chosen  to  record,  and  ill 
f  lit  dll'  process  will  often  be  infinite  in  duration.  Using  only  that  data  we  hav„  recorded  is  equA  - 
alcnt  to  clipping  the  actual  process  at  arbitrary  end  points  in  time.  In  the  frequency  domain  this 
distorts  the  frequency  components  (DFT  coefficients)  from  what  they  would  be  if  we  were  to  consider 
i he  process  is  having  infinite  duration  in  time.  To  reduce  this  distortion,  the  data  points  can  lx 
smoothed  so  that  the  clipping  is  not  so  pronounced.  The  drawback,  however,  is  that  some  of  tlx 
statistical  value  of  the  data  is  lost,  since  smoothing  distorts  the  data  in  the  time  domain.  Tims, 
it  is  wise  to  look  at  both  the  unsmoothed  coefficients  as  wMl  as  the  smoothed  coefficients  in  any 
practical  problem.  In  the  program  described  in  this  paper,  the  smoothed  results  may  either  be 
simply  printed  out  (CALL  SMOOTH!)  or  placed  in  the  data  cells  (CALL  SMOOTH). 

Tin  program  i  ivailable  in  the  form  of  subroutines  assembled  in  a  binary  deck.  It  is  tin 
■  ■pi ion  of  the  user  to  plot  ttie  data  and  to  make  his  own  printouts. 


— -  n(s+l  -j) 

2s 


AJimn.ition  of  the  i*DT  4400  computer  is  that  the  index  on  the  array  H  must  run  from  1  to  N 
in  to.ul  “f  from  otn  >:-i.  Therefore,  in  the  actual  program.  B(j)  is  the  DFT  coefficient  of 
•  ijii<  '!i*  \  (j-I).'NJM  instead  of  j/N  a  t. 


B(n).  n=l .  . . . 
XFORMiM.  B) 

INVERSEfM. 

SMOOT1  l(M.  I 

SMOOTHKM. 


SUMMARY  OF  AVAILABLE  SUBROUTINES 

M 

,  N  is  assumed  to  be  a  complex  array  with  N=2  elements. 

M 

computes  the  2  DF'I  coefficients  of  the  series  B(l).  B(2> . B(N)  and  stores 

these  coefficients  in  the  array  B,  replacing  the  original  series. 

M 

B)  computes  the  2  inverse  OFT  coefficients  of  the  series  B(l),  B(2) . B(N)  and 

stores  these  coefficient-  in  the  array  B.  replacing  the  original  series. 

M 

!)  smooths  the  2  OFT  coefficients  located  in  the  array  B  replacing  the  elements  of 
B  with  these  smoothed  coefficients. 

M 

B)  smooths  the  2  DF1  coefficients  located  in  the  array  B  and  prints  the  smoothed 
coefficients.  The  original  coefficients  are  left  undisturbed. 
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SUBROUTINE  SMOOtHlK.B* 

TYP£  COMPLEX  8.  BEG.  XONE,  XTKO 
DIMENSION  8 ( 102%) 

KPOaa2**x 

8£G«H<1) 

XTROa.25*<-8<2) «2.«B ( '  > -H (KPOa) ) 

*PO*M»KPO»-2 

DO  S  Jal.KPOaM 

XONfcaATaO 

XT*Ob.25*(-B< J>  »?.*B( J«1 )-B  < J»?) ) 

8< J)aXONE 
L  CONTINUE 
XONEaX  T*0 

XTROa.25* ( -8 (KPua-i ) »? ,«B ( KPOB ) -HEG) 
H (XPOa-I ) «XOn£ 

8 (KPOb) bATbO 

HFTURN 

END 


S'lMHOUtlNE  SMOOThI IK,*) 

Type  COMPLEX  8.  BS 
DIMENSION  8(102*) 

POINT  199 
POINT  200 
kP0m*2**K 

8Sa.25* (-8 (2)  *2«*R ( 1 ) -8 (KPOa)  I 
A«L«CaB5(8S) 

JL«n 

AFL«1. 

POINT  201.  JL.  AFL.  A8L.  BS 

KPOaMaKPOa-1 

no  30  J»2«  KPORM 

JL«J-1 

HS«.2S*(-B( JD ♦2.«B(J)-P( J«1 )  ) 

AEL»CABS(HS) 

afl«ael/abl 

POINT  201.  JL.  »FL.  AEl.  BS 
30  CONTINUE 

HS«  . 25* ( “8 ( KPuaM) ♦2,«8(KP0*)-fl(l) ) 

AEL«CAHS(8S) 

AFL* A£L/ A8L 

poInT  201.  KPO«m,  aFL«  AEL,  BS 

190  FOOMAT  (1M1,*  SmOOTHEO  FOURIER  COEFF IC IENTS* »/// » •  REAL  COEE  ■  « 
1 *otAL  Part  OF  CUEF*,/.*  IMAO  COEF  ■  ImaO  part  of  coef*,/. 

?•  AMS  VALUE  a  A«S  VALUE  OF  COEF*./,*  ADJ  COEF  •  A8S  VALUE  • 
3*niVIDE0  BY  ABS  VALUE  OF  COEF  AT  ZERO*.//) 

200  format  (1X,*FR£«J*,3X,*aOJ  C0EF*.3X»*A8S  VALUE*. 3X» 

1  *REAL  COEF*. JX,*lMAG  COEF*,/) 

30)  format  (1X.U,ja,FB.5.3X,E9.3«3X.CCE9.2,Ei2.2)> 

HFTURN 

END 
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APPENDIX  C 


AITI-NOIX  c: 


AN  EXAMPLE 

In  wili-r  !u  show  him  th,  program  works,  we  shall  conipole  llie  l)is,  relc  Fourier  Transform  of 
mi)=siii(-’t  ft)=sin(2ir  I  4>. 


I  ,  i 

N=b4  data  points 
T=t>4 
so  .11  =  I 

k  t  ourier  Series  representation  of  this  function  is  a  pair  of  spikes  at  n=4-lh  and  11--I6,  correspondm 
10  a  sine  wave  of  frequency  f '  -  y  ■  ~  =  -i.  Sampling  x  every  At  —  1  seconds  gives 

\(n)=sin(2  trn/4). 

I  he  modulus  of  each  of  the  l)FT  coefficients  of  x(n)  is  plotted  in  the  accompanying  gTaph.  Note  1 11.1t 
1  In-re  is  indeed  a  peak  at  n=lh.  hut  as  warned  in  appendix  A  there  is  also  a  peak  a.  n=64-l6  4X. 

I  hi--  shows  that  we  can  only  use  the  first  N/2  coefficients  when  trying  to  detect  periodicities.  When 
using  tin.  DFT  as  a  1  ranslorm  in  its  own  right,  'his  restriction  does  not  necessarily  hold. 


OIOS***  TeiT 

Tyoc  CO^uex  A.  CmPl< 

'!|**e'SlON  M10C).  fi(lOO),  X(lCtl) 

r. 

.  *  «eST  bfc  IN  COi-HtEx  ►  0*M  F0»  JSE  IN  blrtPOJTlNh  xrUK* 

TJ  S3  Jilt  04 
«(J)t(J'll*i. 

■jp  »(  )  «CvP.x  (  S  1  \f  (  2  ,  *3  . 14i%V?ob«(  j.l ) /4  ,  ) .  ,  ) 

C  ':i,  slNCc  2«»6«64 

CAI_*  X f  0  4 *  (  6  i  *  I 

(,  »4S»v  4  VC"  CON  T  A  INS  TMfc  DM  COFFHCIEMS 

l' 

r  -s  C1,|,L,c  *ND  PUiT  f H fc  AbbeLUTfe  VALUE  0r  ThF  DM  CO*- ►  H  C  i  ENTs 

"C  4J  J i 1.6 4 

ft1  4(a>«C*hj(*(Jl) 

-1L.  PLC'lrOt x,bi64, - 1 6 , I \Ufc X , 5 , ^COEFFICIENTS, 12- 

1  It*  COEFf  1CIEMS»20 .0  ) 


Kofi’ rollers:  (a)  Cochrm  at  al.  .  "What  is  the  Past  pourin'  I  r  insfoim'"  IHFP  Proceedings, 
pp.  1  004  -1071  Oct  1967 

(h)  Cooky.  Lewis,  and  Welch,  "Application  of  tin  I  act  Pouricr  Transform  to 

Compnt.it ion  of  Fourier  Integrals,  Fourier  Series,  and  Convolution  Integrals,  " 
llvPK  Trans  u  I  ions  on  Audio  and  Fleet  roacoi  1st  ics .  po.  79-84  Jim  1967 

(cl  111  igiian  and  \h.rrow,  ” The  Past  Fourier  I  rnnsforu,  "  IKHK  Spi  ctruiit,  pp.  63-70 
f\  c  i  9o 


Note:  The  June  J967  is  u  of  the  IPF.K  'I  r.insaction  on  Audio  and  h leel roaconsi  ics  is  devoted  to 
applications  of  the  I  a  -t  Fourier  Transform  and  f)i.ere6  huirier  transform. 
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This  research  contribution  describes  a  computer  program  (CNA  Number  76-67)  which 
determines  the  Discrete  Fourier  Transform  of  a  set  of  data,  using  a  recently  developed 
technique  known  as  the  Fast  Fourier  Transform.  The  relation  between  Discrete  Fourier 
Transforms  and  Fourier  Series  when  the  data  is  periodic  is  also  shown. 
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